clc
clear
mm=zeros(20,1);
[D1,S1 ]=xlsread('附件1.csv','附件1','A2:D2227'); %读取基准面主索节点坐标
[D2,S2]=xlsread('附件2.csv','附件2','A2:G2227'); %读取促动器上下端点基准态时端点坐标
[D3,S3]=xlsread('附件3.csv','附件3','A2:C4301'); %读取反射面板主索节点编号
plot3(D1(:,1),D1(:,2),D1(:,3),'b.')  %基准面主索节点坐标
hold on
ffffff=ones(20,1);
for i=1:2226
      plot3([D2(i,1) D2(i,4) ],[D2(i,2) D2(i,5)],[D2(i,3) D2(i,6)],'g') %绘制促动器上下端点连线
      plot3([D2(i,4) D1(i,1) ],[D2(i,5) D1(i,2)],[D2(i,6) D1(i,3)],'y')  %促动器上端点到基准面主索节点
end
XX=mean(D1(:,1));
YY=mean(D1(:,2));
ZZ=mean(D1(:,3));
X2=mean(D2(:,1));
Y2=mean(D2(:,2));
Z2=mean(D2(:,3));
X3=mean(D2(:,4));
Y3=mean(D2(:,5));
Z3=mean(D2(:,6));
plot3([0 XX],[0 YY],[0 ZZ],'r');
inter=10;
for i=1:10
    inter=inter+1;
end
ans=zeros(600,1);%作为备选答案集合
cnt=0;
for g=0.3:0.001:0.6
    f=(139.98+0.33+g)*2;
    c=300+0.4;
    theta=0:0.01:2*pi;
    r=repmat(0:0.5:150,length(theta),1); 
    theta=repmat(theta',1,size(r,2)); 
    x=r.*cos(theta);
    y=r.*sin(theta);
    z=(x.^2+y.^2)/(2*f)-c;
    axis equal
    M=D1(:,1)-D2(:,1); 
    N=D1(:,2)-D2(:,2);
    P=D1(:,3)-D2(:,3);
    k=0;
    DD=zeros(size(D1)); 
    for i=1:length(D1(:,1))
        if (D1(i,1).^2+D1(i,2).^2).^0.5<=150
            k=k+1;
        x0=D2(i,1);
        y0=D2(i,2);
        z0=D2(i,3);
        m=M(i); 
        n=N(i);
        p=P(i);
        if m^2+n^2>0
            t=qiugen(x0,y0,z0,m,n,p,f,c);
            zz(k)=min(z0+p*t);
            t=(zz(k)-z0)/p;
            xx(k)=x0+m*t;
            yy(k)=y0+n*t;
            DD(i,:)=[xx(k) yy(k)  zz(k)];
        else
            xx(k)=x0;
            yy(k)=y0;
            zz(k)=(x0.^2+y0.^2)/(2*f)-c;
            DD(i,:)=[xx(k) yy(k)  zz(k)];
        end
        else
            DD(i,:)=[inf inf inf];
        end
    end
    for i=1:length(D1)
        distance(i)=(sum((DD(i,:)-D1(i,:)).^2))^0.5;
    end
    ll=0;
    distance(find(isinf(distance)))=0; 
    flag=0;
    for kkk=distance
        if(kkk>0.6)
            flag=1;
        end
    end
    if(eq(flag,0))
        cnt=cnt+1;
        ans(cnt)=f;
    end
end
ans=ans./2;%可行的答案集合


%下面进行仿真验证
f=(139.98+0.33+0.4)*2;
c=300.4;
theta=0:0.01:2*pi;  %抛物线极坐标角度范围
r=repmat(0:0.5:150,length(theta),1);   %抛物线极坐标半径范围  定义为矩阵形式
theta=repmat(theta',1,size(r,2));   %抛物线极坐标角度范围 重新调整角度为矩阵形式
x=r.*cos(theta);  %极坐标转化为直坐标
y=r.*sin(theta);   %极坐标转化为直坐标
z=(x.^2+y.^2)/(2*f)-c;  %理想抛物线方程
surf(x,y,z)
shading interp
view([1 0 0])
axis equal
M=D1(:,1)-D2(:,1);  %MNP为空间直线标准方程中间参数
N=D1(:,2)-D2(:,2);
P=D1(:,3)-D2(:,3);
k=0; %口径小于等于300米的范围内主索节点计数
DD=zeros(size(D1)); %初始化记录对应第i个主索节点变化后的位置坐标
for i=1:length(D1(:,1))
    if (D1(i,1).^2+D1(i,2).^2).^0.5<=150  %满足口径小于等于300米的范围要求
        k=k+1;
    x0=D2(i,1); %促动器下端点
    y0=D2(i,2);
    z0=D2(i,3);
    m=M(i);  %以下为直线方程xyz比例系数
    n=N(i);
    p=P(i);
    if m^2+n^2>0   %求出促动器所在直线与理想抛物面的交点
       t=qiugen(x0,y0,z0,m,n,p,f,c);
       zz(k)=min(z0+p*t); %筛选出小于0的Z坐标交点
       t=(zz(k)-z0)/p;
       xx(k)=x0+m*t;
       yy(k)=y0+n*t;
       DD(i,:)=[xx(k) yy(k)  zz(k)];%记录对应第i个主索节点变化后的位置
    else   %如果促动器所在直线与Z轴平行
       xx(k)=x0;
       yy(k)=y0;
        zz(k)=(x0.^2+y0.^2)/(2*f)-c;
      DD(i,:)=[xx(k) yy(k)  zz(k)];%记录对应第i个主索节点变化后的位置
    end
    else
    DD(i,:)=[inf inf inf]; %不在300米口径范围内的则记录为无效值
   end
end
plot3(xx,yy,zz,'r.') %绘制出所有促动器所在直线与理想抛物面的交点
alpha(0.5)


F=[0  0  f/2-c];%理想的抛物线曲面坐标顶点
for i=1:length(DD(:,1)) %依次计算第i个主索节点变化后与方圆1米小圆面的交点

    %P点的Z向坐标值
    PZ=300.4*0.466-300.4; %负值
    %主索节点与焦点方向向量
    mx=DD(i,1);
    ny=DD(i,2);
    pz=DD(i,3)-(f/2-c);
    t=(PZ-(f/2-c))/pz; %根据P点纵坐标求得直线参数方程中间变量
    PX(i)=mx*t;       %得到第i个主索节点变化后与方圆1米小圆面的交点
    PY(i)=ny*t;
    plot3([DD(i,1)  0],[DD(i,2)  0],[DD(i,3)  f/2-c],'b')  %绘制主索节点与焦点连线
end
figure(2) %绘制绘制主索节点与焦点连线与小圆面交点分布图
hold on
for  i=1:length(DD(:,1))
    plot3(PX(i),PY(i),PZ,'k.')  
end
theta=0:0.01:2*pi;  %P点小圆面极坐标角度范围
r=repmat(0:0.05:1,length(theta),1);   %小圆面极坐标半径范围  定义为矩阵形式
theta=repmat(theta',1,size(r,2));   %小圆面v极坐标角度范围 重新调整角度为矩阵形式
x=r.*cos(theta);  %极坐标转化为直坐标
y=r.*sin(theta);   %极坐标转化为直坐标
z=PZ*ones(size(theta));  %纵坐标
surf(x,y,z)
shading interp
alpha(0.5)
view([0  0  1])
            










clc
clear
[D1 S1 ]=xlsread('附件1.csv','附件1','A2:D2227'); %读取基准面主索节点坐标
[D2 S2]=xlsread('附件2.csv','附件2','A2:G2227'); %读取促动器上下端点基准态时端点坐标
[D3 S3]=xlsread('附件3.csv','附件3','A2:C4301'); %读取反射面板主索节点编号
load LL    %导入所有主索节点连接关系矩阵
R=300.4;%基圆半径 可假定为300或300.4两种情况
%绘制理想曲面
tx=279.0:0.05:282;
ty=299.5:0.05:302;   
RD=0:10:130;  %阶梯口径范围
tic
for kk= 1:length(RD)-2
    kk
    clear ra  RLL NLL  rc d1 d2  
 for ii=1:length(tx)
      for jj=1:length(ty)
f=tx(ii);
c=ty(jj);%这里的增量与焦半径的后一个增量相同 %P点向上的部分源于实际数据  过原点抛物线方程的向下平移量  可以调整的量
theta=0:0.01:2*pi;  %抛物线极坐标角度范围
r=repmat(0:0.5:150,length(theta),1);   %抛物线极坐标半径范围  定义为矩阵形式
theta=repmat(theta',1,size(r,2));   %抛物线极坐标角度范围 重新调整角度为矩阵形式
x=r.*cos(theta);  %极坐标转化为直坐标
y=r.*sin(theta);   %极坐标转化为直坐标
z=(x.^2+y.^2)/(2*f)-c;  %理想抛物线方程
%%%%%%%%%%%%%%%%%%求解促动器上下端点连线与理想曲面交点
M=D1(:,1)-D2(:,1);  %MNP为空间直线标准方程中间参数
N=D1(:,2)-D2(:,2);
P=D1(:,3)-D2(:,3);
k=0; %口径小于等于300米的范围内主索节点计数
DD=zeros(size(D1)); %初始化记录对应第i个主索节点变化后的位置坐标
rc=[];
for i=1:length(D1(:,1))
    if (D1(i,1).^2+D1(i,2).^2).^0.5>=RD(kk)&(D1(i,1).^2+D1(i,2).^2).^0.5<=RD(kk+2)  %满足口径范围要求
        rc=[rc i];  %记录在300米口径内的节点序号
        k=k+1;
    x0=D2(i,1); %促动器下端点
    y0=D2(i,2);
    z0=D2(i,3);
    m=M(i);  %以下为直线方程xyz比例系数
    n=N(i);
    p=P(i);
    if m^2+n^2>0   %求出促动器所在直线与理想抛物面的交点
       t=qiugen(x0,y0,z0,m,n,p,f,c);
       zz(k)=min(z0+p*t); %筛选出小于0的Z坐标交点
       t=(zz(k)-z0)/p;
       xx(k)=x0+m*t;
       yy(k)=y0+n*t;
       DD(i,:)=[xx(k) yy(k)  zz(k)];%记录对应第i个主索节点变化后的位置
    else   %如果促动器所在直线与Z轴平行
       xx(k)=x0;
       yy(k)=y0;
        zz(k)=(x0.^2+y0.^2)/(2*f)-c;
      DD(i,:)=[xx(k) yy(k)  zz(k)];%记录对应第i个主索节点变化后的位置
    end
    else
    DD(i,:)=[inf inf inf]; %不在300米口径范围内的则记录为无效值
   end
end



for i=1:length(D1)
    distance(i)=(sum((DD(i,:)-D1(i,:)).^2))^0.5; %直接计算出变化前后的距离
end
distance(find(isinf(distance)))=0;  %将无效距离替换为0
NN(ii,jj,kk)=length(find(distance>0.6));  %显示出不满足距离限制的主索节点个数


%遍历判断口径范围内的反射线是否在1米小圆面以内
%%%%%%%%%%%%%%%%%%%%%主索节点射线与以P为中心点方圆1米小圆面的交点仿真验证
%定义理想抛物曲面的焦点坐标为F
F=[0  0  f/2-c];
pn=0;
Rxy=[];
for i=1:length(rc) %依次计算第i个主索节点变化后与方圆1米小圆面的交点
    %根据P点的Z向坐标值求变化后的第i个主索节点与焦点连线直线方程
    %P点的Z向坐标值
    PZ=R*0.466-R; %负值
    %主索节点与焦点方向向量
    mx=DD(rc(i),1);
    ny=DD(rc(i),2);
    pz=DD(rc(i),3)-(f/2-c);
    t=(PZ-(f/2-c))/pz; %根据P点纵坐标求得直线参数方程中间变量
    PX(i)=mx*t;       %得到第i个主索节点变化后与方圆1米小圆面的交点
    PY(i)=ny*t;
    %plot3([DD(i,1)  0],[DD(i,2)  0],[DD(i,3)  f/2-c],'b')  %绘制主索节点与焦点连线
    if (PX(i)^2+PY(i)^2).^0.5>0.5  %判断主索节点变化后是否在方圆1米小圆面内
        pn=pn+1;
    end
      Rxy=[Rxy (PX(i)^2+PY(i)^2).^0.5]; %记录在方圆1米小圆面内交点的半径
end
 MinRxy(ii,jj,kk)=max(Rxy);  %记录对应本组参数下的最大半径
PN(ii,jj,kk)=pn;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成口径内连接距离变化量分析
LF=zeros(size(LL));  %建立初始判断矩阵
for i=1:length(LL)    %对应所有行
    for j=1:size(LL,2)   %对应两列
          for k=1:length(rc)     %对应300米口径内的第k个节点
               if  LL(i,j)==rc(k)     %如果第k个节点与属于矩阵LL（i，j）相同，则LF取1
                   LF(i,j)=1;
               end
          end
    end
end
lf=sum(LF') ;  %行求和为2则表示该连接存在
RLL=LL(find(lf==2),:); %300米口径内的所有连接关系矩阵
% n=length(RLL) ; %对有效连接进行计数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%根据主索节点调整前后的坐标判断长度变化是否超出了0.07%
%计算RLL连接集合中主索节点调整前D1坐标连接长度与调整后坐标D2连接长度的变化量
for i=1:length(RLL)
     %计算调整前第i个连接的长度
     d1(i)=((D1(RLL(i,1),1)-D1(RLL(i,2),1))^2+(D1(RLL(i,1),2)-D1(RLL(i,2),2))^2+(D1(RLL(i,1),3)-D1(RLL(i,2),3))^2)^0.5;
     %计算调整后第i个连接的长度
     d2(i)=((DD(RLL(i,1),1)-DD(RLL(i,2),1))^2+(DD(RLL(i,1),2)-DD(RLL(i,2),2))^2+(DD(RLL(i,1),3)-DD(RLL(i,2),3))^2)^0.5;
     %计算调整前后变化率
     ra(i)=abs(d2(i)-d1(i))/d1(i);
end
nl(ii,jj,kk)=length(find(ra>0.0007));  %不满足连接长度变化要求的数量
NLL=RLL(find(ra>0.0007),:); %不满足连接长度变化的连接集合
    end
end
%对应3个判断条件  
tmin1=min(min(NN(:,:,kk)))   %口径范围内主索节点满足调节范围判断矩阵
minNN(kk)=tmin1(1);
tmin2=min(min(nl(:,:,kk)))         %口径范围内主索节点连接长度变化判断矩阵
minnl(kk)=tmin2(1);
tmin3=min(min(PN(:,:,kk)))      %口径范围内主索节点穿过小圆面判断矩阵
minPN(kk)=tmin3(1);
end
toc
save result3 NN PN nl  tx ty RD MinRxy



clc      %抛物面方程为x^2+y^2=2p(z-c)时
clear  %txo中保存了分段理想抛物面参数P值； tyo中保存了分段理想抛物面参数C值
%DD中保存了需要调整的主索节点坐标 ,行序号与D1一致
load result3
[D1,S1 ]=xlsread('附件1.csv','附件1','A2:D2227'); %读取基准面主索节点坐标
[D2,S2]=xlsread('附件2.csv','附件2','A2:G2227'); %读取促动器上下端点基准态时端点坐标
[D3,S3]=xlsread('附件3.csv','附件3','A2:C4301'); %读取反射面板主索节点编号

for k=1:length(RD)-2
      F(:,:,k)=nl(:,:,k)+NN(:,:,k)+PN(:,:,k); %生成同时满足3大约束条件的判断矩阵
      n(k)=length(find(F(:,:,k)==0));
end
%求出相邻阶梯式口径范围的可行解集合
for k=1:length(RD)-3
      FF(:,:,k)=F(:,:,k)+F(:,:,k+1);
      nn(k)=length(find(FF(:,:,k)==0));
end
kkkk=ones(30,1);

%查找出一组合理的结果 依次查找前一段区域与后一段区域的可行解交集
for k=1:length(RD)-3
     MinRxy(find(FF(:,:,k)~=0))=100; %把不满足条件的设置为不合理的较大数值
     minr(k)=min(min(MinRxy(:,:,k)))
     [c r]=find(MinRxy(:,:,k)==min(min(MinRxy(:,:,k))));
     txo(k)=tx(c(1));
     tyo(k)=ty(r(1));
end


%绘制分段理想抛物曲面
%先绘制小圆面
figure(1)
    %P点的Z向坐标值
    R=300.4;
    PZ=R*0.466-R; %负值
theta=0:0.01:2*pi;  %P点小圆面极坐标角度范围
r=repmat(0:0.05:1,length(theta),1);   %小圆面极坐标半径范围  定义为矩阵形式
theta=repmat(theta',1,size(r,2));   %小圆面v极坐标角度范围 重新调整角度为矩阵形式
x=r.*cos(theta);  %极坐标转化为直坐标
y=r.*sin(theta);   %极坐标转化为直坐标
z=PZ*ones(size(theta));  %小圆面纵坐标
surf(x,y,z)
shading interp
alpha(0.5)
view([0  0  1])
DD=zeros(size(D1)); %初始化记录对应第i个主索节点变化后的位置坐标
for kk=1:length(RD)-3
f=txo(kk);
 c=tyo(kk);%这里的增量与焦半径的后一个增量相同 %P点向上的部分源于实际数据  过原点抛物线方程的向下平移量  可以调整的量
figure(2)
hold on
 theta=0:0.01:2*pi;  %抛物线极坐标角度范围
r=repmat(RD(kk):0.5:RD(kk+1),length(theta),1);   %抛物线极坐标半径范围  定义为矩阵形式
theta=repmat(theta',1,size(r,2));   %抛物线极坐标角度范围 重新调整角度为矩阵形式
x=r.*cos(theta);  %极坐标转化为直坐标
y=r.*sin(theta);   %极坐标转化为直坐标
z=(x.^2+y.^2)/(2*f)-c;  %理想抛物线方程
surf(x,y,z)
shading interp
view([1 0 0])
%%%%%%%%%%%%%%%%%%求解促动器上下端点连线与理想曲面交点
M=D1(:,1)-D2(:,1);  %MNP为空间直线标准方程中间参数
N=D1(:,2)-D2(:,2);
P=D1(:,3)-D2(:,3);
k=0; %口径小于等于300米的范围内主索节点计数
rc=[];
for i=1:length(D1(:,1))
    if (D1(i,1).^2+D1(i,2).^2).^0.5<=RD(kk+1)&(D1(i,1).^2+D1(i,2).^2).^0.5>=RD(kk) %满足口径小于等于300米的范围要求
        rc=[rc i];  %记录在300米口径内的节点序号
        k=k+1;
    x0=D2(i,1); %促动器下端点
    y0=D2(i,2);
    z0=D2(i,3);
    m=M(i);  %以下为直线方程xyz比例系数
    n=N(i);
    p=P(i);
    if m^2+n^2>0   %求出促动器所在直线与理想抛物面的交点
       t=qiugen(x0,y0,z0,m,n,p,f,c);
       zz(k)=min(z0+p*t); %筛选出小于0的Z坐标交点
       t=(zz(k)-z0)/p;
       xx(k)=x0+m*t;
       yy(k)=y0+n*t;
       DD(i,:)=[xx(k) yy(k)  zz(k)];%记录对应第i个主索节点变化后的位置
    else   %如果促动器所在直线与Z轴平行
       xx(k)=x0;
       yy(k)=y0;
        zz(k)=(x0.^2+y0.^2)/(2*f)-c;
      DD(i,:)=[xx(k) yy(k)  zz(k)];%记录对应第i个主索节点变化后的位置
    end
   end
end
plot3(xx,yy,zz,'r.') %绘制出所有促动器所在直线与理想抛物面的交点
alpha(0.5)


% %%%%%%%%%%计算经过主索节点射线后与方圆1米小圆面的交点
ttt=0;
for i=1:100
    ttt=ttt+1;
end
kkk=0;
for i=1:100
    kkk=kkk+1;
end
ggg=0;
for i=1:100
    ggg=ggg+1;
end

ldjald=1;
% figure(2)
F=[0  0  f/2-c];
for i=1:length(DD(:,1)) %依次计算第i个主索节点变化后与方圆1米小圆面的交点
    mx=DD(i,1);
    ny=DD(i,2);
    pz=DD(i,3)-(f/2-c);
    t=(PZ-(f/2-c))/pz; %根据P点纵坐标求得直线参数方程中间变量
    PX(i)=mx*t;       %得到第i个主索节点变化后与方圆1米小圆面的交点
    PY(i)=ny*t;
end
figure(1) %绘制绘制主索节点与焦点连线与小圆面交点分布图
hold on
for  i=1:length(DD(:,1))
    plot3(PX(i),PY(i),PZ,'k.')  
end
end


%%%反射面计数

ggg=zeros(60,1);
SD=sum(DD');
rc=find(SD~=0); %节点序号
jadjal=zeros(30,1)
ahfka=0;
for i =jadjal
    ahfka=ahfka+1;
end

DS=zeros(size(S3)); %初始化判断矩阵  用于判断反光板是否在300米口径内
for i=1:size(S3,1)              %按行检索S3
       for  j=1:size(S3,2)      %按列检索S3
           for k=1:length(rc)  %依次比对S3{i,j}是否属于rc
                  if strcmp(S3{i,j},S1{rc(k)})    %如果S3{i,j}属于rc，则DS(i,j)取1
                    DS(i,j)=1;  
                    break
                  end
           end
       end
end
tf=sum(DS') ;  %判断三角反射板三个顶点是否都在300米口径内
ns=length(find(tf==3));  %对有效三角反射板进行计数
kuan=0;
if(ns>10000)
    kuan=100000;
end

